function [h0Hat,hxHat,sigmaHat,scalingOpt] = runInduceStatVAR(hx,factors,VarFactors,Cov10_u,Cov01_u)

T       = size(factors,2);
modulus = 100;
if any(modulus>=1)
    scaling = 1;
    while any(modulus>=1)
        scaling       = scaling*0.99;
        hxTmp         = scaling*hx;
        eigValues     = eig(hxTmp);
        realEigValues = real(eigValues);
        imagEigValues = imag(eigValues);
        modulus       = sqrt(realEigValues.^2 + imagEigValues.^2);
    end
    % Moment based scaling to induce stationarity
    meanFactors   = mean(factors,2);
    input.varXData= (factors-repmat(meanFactors,1,T))*(factors-repmat(meanFactors,1,T))'/(T-1)-mean(VarFactors,3);  
    input.Var_u   = VarFactors;
    input.Cov10_u = Cov10_u;
    input.Cov01_u = Cov01_u;
    input.factors = factors;
    maxIter       = 10000;
    tolFun        = 1e-8;
    options       = optimset('Display','off','MaxIter',maxIter,'MaxFunEvals',10000,'TolFun',tolFun,'TolX',1e-8);
    input.hx      = hx;
    % Nelder-Mead optimization
    scalingOpt = fminsearch(@induceStatVAR,scaling,options,input);
    if 1 == 0
        % Newton Raphson optimization - implemented in the fortran version
        diffQ    = 1;
        idxIter  = 1;
        epsValue = 1e-6;
        while idxIter <= maxIter && abs(diffQ) > tolFun
            % The gradient
            Q_p_eps  = induceStatVAR(scaling+epsValue,input);
            Q_m_eps  = induceStatVAR(scaling-epsValue,input);
            grad     = (Q_p_eps - Q_m_eps)/(2*epsValue);
            % The Hessian
            Q        = induceStatVAR(scaling,input);
            Q_pp_eps = induceStatVAR(scaling+2*epsValue,input);
            Q_mm_eps = induceStatVAR(scaling-2*epsValue,input);
            hess     = ((Q_pp_eps-Q)-(Q-Q_mm_eps))/(4*epsValue^2);
            hess     = abs(hess);   %ensure that hess is positive!
            % Step
            delta      = - grad/hess;
            scalingNew = scaling + delta;
            
            % Updating
            QNew       = induceStatVAR(scalingNew,input);
            diffQ      = QNew-Q;
            idx        = 0;
            while diffQ > 0 && idx < 100
                % We use a scaling of delta 
                idx        = idx + 1;
                scalingNew = scaling + 0.80^idx*delta;
                QNew       = induceStatVAR(scalingNew,input);
                diffQ      = QNew-Q;
            end
            scaling    = scalingNew;
            idxIter    = idxIter + 1;
            %[QNew, idx, idxIter, scalingNew]
        end
        scalingNewtonRaphson = scaling;
    end
    [Qopt,h0Hat,hxHat,sigmaHat] = induceStatVAR(scalingOpt,input);
    disp(['Objective function when inducing stationarity of VAR model : Qopt = ',num2str(Qopt)]);
end

end

